/* ------------------------------------------------------------------
		   Support code for problem set #5

			     SOURCE CODE
------------------------------------------------------------------- */

#include "support.hh"
#include <stdlib.h>
#include <sys/types.h>
#include <fstream.h>
#include <math.h>

ifstream parameter_file("parameters");

void read_data(string filename, matrix& input, matrix& output)
{
  ifstream data(filename);

  int patterns, inputs, outputs;

  data.scan("# patterns = %d", &patterns);
  data.scan("# inputs = %d", &inputs);
  data.scan("# outputs = %d", &outputs);

  input.resize(patterns, inputs);
  output.resize(patterns, outputs);

  for (int i = 0; i < patterns; i++) {

    for (int j = 0; j < inputs; j++)
      data >> input[i][j];

    for (j = 0; j < outputs; j++)
      data >> output[i][j];
  }
}

int seeded = 0;

double random_number()
{
  if (!seeded) {
    srandom(getpid());
    seeded = 1;
  }

  return ((double) (random() % 100000)) / 50000.0 - 1.0;
}

char fscan_until(istream &is, string& str)
{
  if (str.len() == 0 ) return 0;

  char lastc = 0, c;
  int i = 0;

  do {
    is.get(c);
    if (lastc == '\n' && c == '#') { 
      i = 0; 
      fscan_until(is, "\n");
      lastc = '\n';
    }
    else {
      if (c == str[i])
	i++;
      else if (c == str[0])
	i = 1;
      else 
	i = 0;

      lastc = c;
    }
  }
  while (str[i] != '\0' && !is.eof());

  return is.eof();
}

void read_parameter(string& parameter_name, double& parameter, 
		    double default_value)
{
  if (parameter_file.bad()) {
    cerr << "ERROR: Attempt to read nonexistent parameter file\n";
  }

  parameter_file.seekg(0);
  fscan_until(parameter_file, parameter_name);

  parameter_file >> parameter;

  if (parameter_file.fail())
    parameter = default_value;
}

void read_parameter(string& parameter_name, int& parameter, 
		    int default_value)
{
  if (parameter_file.bad()) {
    cerr << "ERROR: Attempt to read nonexistent parameter file\n";
  }

  parameter_file.seekg(0);
  fscan_until(parameter_file, parameter_name);

  parameter_file >> parameter;

  if (parameter_file.fail())
    parameter = default_value;
}

vector augment(vector& v) return w(v);
{
  w.upsize(v.dimension() + 1);
  w[v.dimension()] = 1;
}

class SVD_matrix {
/* ===================================================================
			P   U   B   L   I   C
=================================================================== */
public:

  SVD_matrix(const matrix& m, double tol = 1E-12);

  friend SVD_matrix SVD(const matrix& m);
  friend SVD_matrix SVD(const matrix& m, double tol);  

  vector solve(vector b) const;
  vector solve(vector b, vector& residual) const;
  matrix pseudoinverse(void) const;

  vector singular_values(void) { return w; }
  matrix range_basis(void) {return u; }
  matrix kernel_basis(void) { return v; }

  void set_tolerance(double new_tol) { tolerance = new_tol; }
  
/* ===================================================================
		      P  R  O  T  E  C  T  E  D
=================================================================== */
protected:

  double tolerance;                // zero for small singular values
  int n_rows;
  int n_cols;
  matrix u;                 // Columns of U give range basis
  vector w;                 // Singular values
  matrix v;                 // Columns of V give kernel basis

  double pythag(double a, double b);
  void SVD_decomp(const matrix& mat);
};

inline SVD_matrix SVD(const matrix& m) return svd(m);
{
}

inline SVD_matrix SVD(const matrix& m, double tol)
     return svd(m, tol);
{
}

matrix pseudoinverse(matrix& M)
{
  SVD_matrix svd(M);
  return svd.pseudoinverse();
}

/******************************************************************/
/*                     Numerical Recipes Macros                   */
/******************************************************************/
inline float SQR(float a) {
  return (a == 0.0 ? 0.0 : a * a);
}
inline double SQR(double a) {
  return (a == 0.0 ? 0.0 : a * a);
}

inline double MAX(double a, double b) {
  return (a > b ? a : b);
}
inline double MIN(double a, double b) {
  return (a < b ? a : b);
}

inline float MAX(float a, float b) {
  return (a > b ? a : b);
}
inline float MIN(float a, float b) {
  return (a < b ? a : b);
}

inline int MAX(int a, int b) {
  return (a > b ? a : b);
}
inline int MIN(int a, int b) {
  return (a < b ? a : b);
}

inline long MAX(long a, long b) {
  return (a > b ? a : b);
}
inline long MIN(long a, long b) {
  return (a < b ? a : b);
}

inline double SIGN(double a, double b) {
  return  (b >= 0.0 ? fabs(a) : -fabs(a));
}

// return a square matrix whose diagonal elements are given by v
matrix diag(const vector& v)
     return m(v.dimension(), v.dimension());
{
  for (register int i=0; i<v.dimension(); i++)
    m[i][i] = v[i];
}

vector diag(const matrix& mat)
{
  int r = mat.rows();
  int c = mat.columns();

  if (r != c) {
    cout << "ERROR: Nonsquare matrix in routine diag(matrix)\n";
    exit(0);
  }
  vector diag_v(r);

  for (register int i=0; i<r; i++)
    diag_v[i] = mat[i][i];

  return diag_v;
}

matrix diagonal_entries(const matrix& mat)
{
  int r = mat.rows();
  int c = mat.columns();

  if (r != c) {
    cout << "ERROR: Nonsquare matrix in method diagonal_entries\n";
    exit(0);
  }
  matrix diag_mat(r, c);

  for (register int i=0; i<r; i++)
    diag_mat[i][i] = mat[i][i];

  return diag_mat;
}


SVD_matrix::SVD_matrix(const matrix& m, double tol)
: u(m.rows(), m.columns()), w(m.columns()), v(m.columns(), m.columns())
{
  tolerance = tol;
  n_rows = m.rows();
  n_cols = m.columns();
  SVD_decomp(m);
}

vector SVD_matrix::solve(vector b) const return x(n_cols);
{
  // Pseudoinverse is an n_cols by n_rows matrix
  x = pseudoinverse() * b;
}

vector SVD_matrix::solve(vector b, vector& residual) const return x(n_cols);
{
  x = pseudoinverse() * b;
  // residual = mx - b
  residual = u * diag(w) * transpose(v) * x - b;
}

matrix SVD_matrix::pseudoinverse(void) const
     return pseudo(n_cols, n_rows);
{

  // Pseudoinverse is v * w^-1 * u^T

  // Pseudoinvert the diagonal matrix of singular values
  vector w_inv(w);  
  for (register int i=0; i<n_cols; i++) {
    if (w[i]) w_inv[i] = 1.0/w[i];
    else w_inv[i] = 0.0;
  }

  pseudo = v * diag(w_inv) * transpose(u);
}

/* -------------------------------------------------------------------
			  Protected methods
------------------------------------------------------------------- */

/*     (C) Copr. 1986-92 Numerical Recipes Software #.3.   */

double SVD_matrix::pythag(double a, double b)
{
  double absa = fabs(a);
  double absb = fabs(b);

  if (absa > absb) return absa * sqrt(1.0 + SQR(absb/absa));
  else return (absb == 0.0 ? 0.0 : absb * sqrt(1.0 + SQR(absa/absb)));
}

void SVD_matrix::SVD_decomp(const matrix& mat) {

  matrix a(mat);

  int m = n_rows;
  int n = n_cols;

  int flag,i,its,j,jj,k,l,nm;
  double anorm, c, f, g, h, s, scale, x, y, z;

  vector rv1(n);

  g=scale=anorm=0.0;

  // Begin Householder reduction to bidiagonal form
  for (i=0;i<n;i++) {
    l=i+1;
    rv1[i]=scale*g;
    g=s=scale=0.0;
    if (i < m) {
      for (k=i;k<m;k++) scale += fabs(a[k][i]);
      if (scale) {
	for (k=i;k<m;k++) {
	  a[k][i] /= scale;
	  s += a[k][i]*a[k][i];
	}
	f=a[i][i];
	g = -SIGN(sqrt(s),f);
	h=f*g-s;
	a[i][i]=f-g;
	for (j=l;j<n;j++) {
	  for (s=0.0,k=i;k<m;k++) s += a[k][i]*a[k][j];
	  f=s/h;
	  for (k=i;k<m;k++) a[k][j] += f*a[k][i];
	}
	for (k=i;k<m;k++) a[k][i] *= scale;
      }
    }
    w[i]=scale * g;
    g=s=scale=0.0;
    if (i < m && i != n-1) { 
      for (k=l;k<n;k++) scale += fabs(a[i][k]);
      if (scale) {
	for (k=l;k<n;k++) {
	  a[i][k] /= scale;
	  s += a[i][k]*a[i][k];
	}
	f=a[i][l];
	g = -SIGN(sqrt(s),f);
	h=f*g-s;
	a[i][l]=f-g;
	for (k=l;k<n;k++) rv1[k]=a[i][k]/h;
	for (j=l;j<m;j++) {
	  for (s=0.0,k=l;k<n;k++) s += a[j][k]*a[i][k];
	  for (k=l;k<n;k++) a[j][k] += s*rv1[k];
	}
	for (k=l;k<n;k++) a[i][k] *= scale;
      }
    }
    anorm=MAX(anorm,(fabs(w[i])+fabs(rv1[i])));
  }  // End Householder section

  // Accumulation of right-hand transformations
  for (i=n-1;i>=0;i--) {
    if (i < n-1) {
      if (g) {
	for (j=l;j<n;j++)     // Double division to avoid underflow
	  v[j][i]=(a[i][j]/a[i][l])/g; 
	for (j=l;j<n;j++) {
	  for (s=0.0,k=l;k<n;k++) s += a[i][k]*v[k][j];
	  for (k=l;k<n;k++) v[k][j] += s*v[k][i];
	}
      }
      for (j=l;j<n;j++) v[i][j]=v[j][i]=0.0;
    }
    v[i][i]=1.0;
    g=rv1[i];
    l=i;
  } // End right-hand transformation accumulation

  // Accumulation of left-hand transformations
  for (i=MIN(m-1,n-1);i>=0;i--) {
    l=i+1;
    g=w[i];
    for (j=l;j<n;j++) a[i][j]=0.0;
    if (g) {
      g=1.0/g;
      for (j=l;j<n;j++) {
	for (s=0.0,k=l;k<m;k++) s += a[k][i]*a[k][j];
	f=(s/a[i][i])*g;
	for (k=i;k<m;k++) a[k][j] += f*a[k][i];
      }
      for (j=i;j<m;j++) a[j][i] *= g;
    } else for (j=i;j<m;j++) a[j][i]=0.0;
    ++a[i][i];
  } // end left-hand transformation accumulation

  // Diagonalization of the bidiagonal form
  for (k=n-1;k>=0;k--) {
    for (its=1;its<=30;its++) {    // Allow 30 iterations
      flag=1;
      for (l=k;l>=0;l--) {         // Test for splitting
	nm=l-1;                    // Note that rv1[0] is always 0.0
	if ((double)(fabs(rv1[l])+anorm) == anorm) {
	  flag=0;
	  break;
	}
	if (nm < 0) cout << "WARNING: nm out of range for a vector\n";
	if ((double)(fabs(w[nm])+anorm) == anorm) break;
      }
      if (flag) {
	c=0.0;                     // Cancellation of rv1[l], if l > 0
	s=1.0;
	for (i=l;i<=k;i++) {
	  f=s*rv1[i];
	  rv1[i]=c*rv1[i];
	  if ((double)(fabs(f)+anorm) == anorm) break;
	  g=w[i];
	  h=pythag(f,g);
	  w[i]=h;
	  h=1.0/h;
	  c=g*h;
	  s = -f*h;
	  for (j=0;j<m;j++) {
	    y=a[j][nm];
	    z=a[j][i];
	    a[j][nm]=y*c+z*s;
	    a[j][i]=z*c-y*s;
	  }
	}
      }
      z=w[k];
      if (l == k) {                 // Convergence
	if (z < 0.0) {              // Singular value is made nonnegative
	  w[k] = -z;
	  for (j=0;j<n;j++) v[j][k] = -v[j][k];
	}
	break;
      }
      if (its == 30) cout << "no convergence in 30 svdcmp iterations\n";;
      x=w[l];                       // Shift bottom 2x2 minor
      nm=k-1;
      y=w[nm];
      g=rv1[nm];
      h=rv1[k];
      f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
      g=pythag(f,1.0);
      f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
      c=s=1.0;                      // Next QR transformation
      for (j=l;j<=nm;j++) {
	i=j+1;
	g=rv1[i];
	y=w[i];
	h=s*g;
	g=c*g;
	z=pythag(f,h);
	rv1[j]=z;
	c=f/z;
	s=h/z;
	f=x*c+g*s;
	g = g*c-x*s;
	h=y*s;
	y *= c;
	for (jj=0;jj<n;jj++) {
	  x=v[jj][j];
	  z=v[jj][i];
	  v[jj][j]=x*c+z*s;
	  v[jj][i]=z*c-x*s;
	}
	z=pythag(f,h);
	w[j]=z;                      // Rotation can be arbitrary if z=0
	if (z) {
	  z=1.0/z;
	  c=f*z;
	  s=h*z;
	}
	f=c*g+s*y;
	x=c*y-s*g;
	for (jj=0;jj<m;jj++) {
	  y=a[jj][j];
	  z=a[jj][i];
	  a[jj][j]=y*c+z*s;
	  a[jj][i]=z*c-y*s;
	}
      }
      rv1[l]=0.0;
      rv1[k]=f;
      w[k]=x;
    }
  }
  u = a;
  // set small singular values to zero
  for (i = 0; i< n; i++)
    if (w[i] < tolerance) w[i] = 0.0;
}